CartPoleを機械学習じゃなくて古典制御でやってみた

#Python

山名

2024.05.17

はじめに

皆さんはCartPoleをご存知でしょうか。CartPoleはOpenAI社が機械学習の勉強やアルゴリズムの開発のために提供している(していた)オープンソースのライブラリに含まれるゲームの1つです。2021年にOpenAIがGymの更新をやめてからは、非営利団体によってforkされた系のGymnasiumが更新されています。

さて、上述のようにCartPoleは機械学習で遊ぶものですが、古典制御でも同様の倒立制御ができるような気がしてきました。早速やってみます。

古典制御にもいくつか種類がありますが、今回は一番有名なPID制御を用いることにします。ただし、伝達関数を定義して解くことはせず、観察しながら試行錯誤でゲインを決定していこうと思います。(シミュレーター上の倒立振子の伝達関数を求め出したらそれだけでブログを5本くらい書けちゃいそうなので)

※本記事内ではGymnasiumのインストール部分やCartPoleについての細かいコードの説明はしません。

※グラフによって軸のスケールが異なるのでご注意ください。

※学生時代に伝達関数の計算が難しすぎて試験前に泣いていました。

乱数

まずは比較用に乱数で制御します。main.pyは以降共通でエントリーポイントになります。

from my_gym import Gym

g = Gym(200)
g.act()
g.output()

import gymnasium as gym
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

class Gym():
    def __init__(self, t):
        self.env = gym.make('CartPole-v1', render_mode='rgb_array')
        self.T = t

        # 初期化
        self.state, self.info = self.env.reset()
        self.state_data = []
        self.render_data = [self.env.render()]

    def act(self):
        for t in range(self.T):
            # 乱数で移動
            action = self.env.action_space.sample()
    
            # 状態を遷移
            next_state, reward, terminated, truncated, _ = self.env.step(action)

            if terminated is True or truncated is True:
                break
            
            # 状態を保存
            self.state_data.append((self.state, action, reward, terminated))
            self.render_data.append(self.env.render())

            # 更新
            self.state = next_state

        # 終了時の状態
        self.state_data.append((self.state, None, None, None))
    
    def output(self):
        def update(t):
            # 状態を取得
            rgb_data = self.render_data[t]

            # 描画
            plt.imshow(rgb_data)
            plt.xticks(ticks=[])
            plt.yticks(ticks=[])

        # gifを作成
        fig = plt.figure(figsize=(9, 7), facecolor='white')
        anime = FuncAnimation(fig=fig, func=update, frames=len(self.render_data), interval=50)
        anime.save('CartPole_random.gif')

結果

ステップ数は19でした。ランダムに動いたらこんなもんですね。中途半端なところで終わっているように見えますが、CartPoleの終了条件の1つに「角度が±12°を超える」があるためです。

P制御

続いてP制御をやってみます。P制御では、目標値と現在値の偏差を0に近づけるように制御していきます。

# p制御
def P(current, target=0):
    if current < target:
        return 0 # left
    return 1 # right

こんな関数を作りました。CartPoleでは、振子の角度が取れるのでこれを0[rad]に近づければ良さそうです。my_gym.pyの19行目を、

# action = self.env.action_space.sample()
# ↓
action = P(self.state[2])

こんなふうに変えます。state[2]には振子の角度が入っています。

結果

ステップ数は56でした。微成長です。オーバーシュートしていて、掃除中に手のひらに立てたホウキを追いかけて走った記憶が思い出されます。

ところで、CartPoleのx軸の移動は「左は0、右は1」の2値で制御しています。ゲインが介在していないのは制御として微妙な気がします。

改良型P制御

幸いなことに、Gymnasiumのコードベースは読める量に収まっています。ライブラリ中の以下を改変してx軸方向の変位を数値として入力できるようにしてみました。コメントの行番号は0.29.1の内容です。

def step(self, action): # 131行目
    # assert self.action_space.contains(        #　ここ
    #     action                                #　ここ
    # ), f"{action!r} ({type(action)}) invalid" #　ここ
    assert self.state is not None, "Call reset before using step method."
    x, x_dot, theta, theta_dot = self.state
    force = action # ここ
    costheta = math.cos(theta)
    sintheta = math.sin(theta)

※GymnasiumはMITライセンスです。

ライブラリ側の変更に合わせてコード側も変更していきます。ゲインに当たる変数:Kpの値を調整しながらできるだけ倒立を維持できるように試行を繰り返します。

def P(current, target=0):
    Kp = 0.125
    # 100を掛けてcartの移動のスケールに合わせる
    return -((target-current)*Kp*100)

何回か試行してKpの値が0.125くらいがちょうどいい感じがしました。このときの各要素の遷移が以下のグラフです。

結果

また、ビジュアライズするとこうなります。

倒立は維持していますが、cartがどんどんずれていることがわかります。これ以上Kpを大きくするとオーバーシュートが起こり、これ以上小さくすると角度を反転させるに足る力を得られません。P制御だけだとこんなもんな気がします。（cartの位置に制約がなければ収束させることもできるかも?）

また、ここら辺から結果が見え始めるのが遅くなってくるのでステップ数を500にして、10ステップごとに描画しています。

PD制御

名前的には次にI制御を足したいところですが、観察した結果を見ると今必要なのはオーバーシュートの抑制であるので、まずはD制御を足します。観察しながらパラメータを決定してもいいんですが、せっかくなのでジーグラニコルスの限界感度法で決定してみます。この手法は、P制御でちょうど振動するときのKpをKu、振動の周期をPuとして下記表の定数をつける方法です。

パラメータをいじって振動するKpを探します。Ku=0.12くらいで若干発散に向かっている気もしますが振動しました。

ここから、

パラメータ	値
Kp	0.6Ku
Pu	200
Ti	0.5Pu
Td	0.125Pu

の各値を導けます。PD制御のパラメータは厳密には出せませんが、以上を踏まえて、

Ku = 0.12
Pu = 200
Kp = Ku * 2 # このくらいの大きさにしないとうまくいかない
def P(current, target=0):
    return -((target-current)*Kp*100)

def D(current_speed, target=0):
    Kd = Kp * Pu * 0.3
    return -Kd * (target - current_speed)

試行錯誤の上このように決定し、current_sppedにself.state[3]を入力して実行します。反応が遅く大きいので、ステップ数は1000回、終了フラグは無視して試行しました。結局ほぼ試行錯誤でパラメータを決定しています。限界感度法は通用しなかったけどそれっぽい気分は味わえるのでOKです。

ちなみに、D制御のパラメータとしてspeedを与えていますが、角速度が角度の時間微分であるためです。

結果

オーバーシュートを抑制して倒立できていますが、迷いなく場外に逃げていきます。これはcartの位置に関しての制御をしていないためと考えられます。そこで、若干悪手な気もしますが中心に留まるようなP制御とD制御を足してみます。currentにはstate[0], current_speedにはstate[1]を代入します。

def Px(current, target=0):
    return -((target-current))

def Dx(current_speed, target=0):
    Kd = Kp * Pu * 0.1
    return -Kd * (target - current_speed)

これで実行すると、

結果

静止画かと思うほど直立しています。とても良さそうです。

PID制御

いよいよPID制御をやっていきます。I制御を足すことで、一定時間正負どちらかに振れ続けるようなときに強い反対向きの力を加えることができます。例によってcartの位置にもI制御を足します。~~正直すでにうまいこといってるので、これで終わりでもいいんですが。~~

ここまでの関数も含めるとこんな感じです。

Ku = 0.12
Pu = 200
Kp = Ku * 2
def P(current, target=0):
    return -((target-current)*Kp*100)

def Px(current, target=0):
    return -((target-current)*1)

def I(t, state_data, target=0):
    Ti = 0.3 * Pu
    sum = 0
    for i in range(t):
        sum += target-state_data[i][0][2]
    Ki = Kp/Ti
    return -Ki * sum * 100

def Ix(t, state_data, target=0):
    Ti = 0.9 * Pu
    sum = 0
    for i in range(t):
        sum += target-state_data[i][0][0]
    Ki = Kp/Ti
    return -Ki * sum

def D(current_speed, target=0):
    Kd = Kp * Pu * 0.3
    return -Kd * (target - current_speed)

def Dx(current_speed, target=0):
    Kd = Kp * Pu * 0.1
    return -Kd * (target - current_speed)

時間での0->tの定積分として総和を用いています。本来は連続データなので意味が変わってきそうですが、listに入ってるなら離散データと言える気もするし大きく食い違っているわけでもないのでいいかなと思っています。何度か試してこのくらいのパラメータがちょうど良さそうでした。気になる結果は...